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A spin dynamics algorithm, combining checkerboard updating and a rotation algorithm based on 
the local second-rank ordering field, is developed for the Lebwohl-Lasher model of liquid crystals. 
The method is shown to conserve energy well and to generate simulation averages which are con- 
sistent with those obtained by Monte Carlo simulation. However, care must be taken to avoid the 
undesirable effects of director rotation, and a method for doing this is proposed. 

PACS numbers: 05.10.-a,61.20. Ja,61.30.Cz,64.70.Md 

Many physical systems may be represented by the highly idealized model of a set of interacting classical spin vectors 
located on a regular, often cubic, lattice lj. This paper considers a classical Hamiltonian of the form 

H=J2 H Sj = £ h(s 3 ■ s k ) (1) 
(j.k) j keNj 

where the interaction energy is given by 

h( Sj ■ s k ) = -J (|( Sj • s k f - §) (2) 

In eqn (Q), the notation k ^ indicates a sum over all nearest neighbour lattice sites j, k, considering each neighbour 
pair only once. Correspondingly J^keN denotes a sum over sites k which constitute the set Nj of nearest neighbours 
of j. Full periodic boundary conditions are assumed. The spins s 3 are three-dimensional vectors of unit length, and 
the pair interaction has the form of a simple function h of the scalar product of the corresponding vectors. The case 
h(sj ■ Sk) — —Jsj ■ Sk is the well-studied classical Heisenberg model, which, for J > 0, exhibits a ferromagnetic ground 
state. Our interest lies in the model defined by eqn J2Jl, originally proposed by Lebwohl and Lasher 0,|j|to represent 
the ordering in nematic liquid crystals. This model has been extensively studied by Monte Carlo simulations |l| in 
the canonical ensemble 0, 0, El • At high temperature T the system is disordered, while for J > the low-temperature 
phase has aligned spins; the definition of an order parameter will be given below. The phase transition is known to 
be weakly first order. Brownian or Langevin dynamics have also been applied to this model [UBS El; here the 
approach of spin dynamics is considered. 

The spin-dynamics equations of motion take the form 

. dH 

S j — ~~~ X S j — & ' j X S j ( o J 

dsj 

the dot denotes time differentiation, and x is the vector product. The local field dH/dsj at spin j defines an 
instantaneous angular velocity of precession, Slj. This dynamics conserves the individual spin lengths, as may be seen 
by considering the time derivative of \sj\ 2 : 

^\s 3 f = 2 Sj ■ Sj = 2(Slj x Sj ) -Sj=0. 



The hamiltonian is also conserved: 

ds 
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Finally, for the class of hamiltonians of eqn QJ, the total magnetization S = Y^j s j is a l so conserved: 

dH 
dsj 



j 3 j fceiVj 



Here h' stands for the derivative of h with respect to its argument, and the right hand side vanishes because every 
pair interaction is included twice, once as (j, k) and once as (k,j), and these cancel because S X S j — Sj X S . 
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An algorithm to simulate the spin dynamics of the Heisenberg model was developed independently by Frank et al. 
[Tl| and Krech et al. The set of all spins s = {sj} is subdivided into two sublattices, A and B, in a checkerboard 
fashion. All spins on sublattice A interact only with their nearest neighbours, which are all on sublattice B, and vice 
versa. Formally, the generator of infinitesimal rotations of the whole set of spins may be decomposed into separate, 
non-commuting, matrices or operators which act on the corresponding sublattices. A time-reversible approximation 
to the full rotation operator acting over a finite timestep may be formally written 

s(At) « B(±At)A(At)B(±At)s + 0{At 3 ) (4) 

The operator A(At) represents the rotation of the A-spins in an external field caused by the neighbouring fixed 
B-spins, during a period At; similarly for B. The algorithm proceeds in an alternating fashion: first solving the 
dynamics of the spins on one sublattice, with the other sublattice spins held fixed; then vice versa, and so on. For 
the detailed justification of this algorithm, see [TlUT^] . 

For the Heisenberg model, this decomposition is especially convenient, since Q.j does not depend on Sj, and the 

motion of spins on a sublattice Sj = Clj x Sj with constant Slj may be solved exactly. Set Clj = \ ftj\, and tlj = Qj/ilj 
in the finite rotation formula |l3j to give 

Sj(At) = (Clj ■ Sj)Sij + sin(fij-Ai) Clj x Sj + cos(fijAi) [Sj - (Clj ■ s^fl^ . (5) 

This represents the practical implementation of one of the A or B sub-steps in eqn . 
For the Lebwohl-Lasher model of eqn J5J , 

ftj = -3 J 2J ( S J ■ s k)sk , 

keN-i 

the motion during one time step is no longer a simple rotation, because although the neighbouring spins Sfc are fixed, 
the moving spin Sj appears on the right. Krech et al. [T^ propose an iterative approach to this problem; here a 
different method is adopted. The equation of motion of each spin is conveniently written 

Sj = il x Sj = -3 J ( s i ■ x s i = S J" ' F - x S J 

keNj 

defining a tensor field at each lattice site due to the neighbouring spins 

Fj = -3J Y ( s fe ® s k - I 1 ) ■ 

A term ^1, where 1 is the unit tensor, is subtracted to make Fj traceless: this has no effect on the equations of motion 
since Sj ■ 1 x Sj = Sj x Sj = 0. 

The above equation is well known in another context: the torque-free time evolution of the angular velocity of a 
rigid body, expressed in the frame of reference of the inertia tensor of the body itself |l3l |. The symmetric tensor 
Fj plays the role of the inertia, but it arises from a different mechanism here. A method for integrating this over a 
timestep has been proposed 14]. It is convenient to resolve all the vectors in the principal axis system of Fj. Denote 
the eigenvalues, in order Fji > Fj2 > Fj%, and the corresponding eigenvectors lj,2j,3j. These are taken to be 
mutually perpendicular and of unit length, and all the following vector and matrix expressions are expressed in this 
frame. It should be noted that the algorithm is independent of the sign ambiguity associated with these eigenvectors. 
Fj becomes diagonal, and the instantaneous angular velocity takes a simple form: 

(Fj X \ (sjiFjA 
Fj = F j2 , => Slj = s j2 F j2 
V F j3 J \Sj3F jS J 

The equations of motion Sj = ilj x Sj become 

Sji = [Fj2 — Fj3) s j2 s j3 ! an d cyclic permutations. 

Consider the component involving Fji, generating rotations about the lj axis (the others are similar): 



Sji - , Sj 2 - -F,jiSji Sj 3 , Sj 3 - FjlSji s j2 
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During the action of this part of the field, the component Sj\ remains constant, and the other two components are 
rotated about the lj direction at an angular velocity = FjiSj±. Suppose the generator of this infinitesimal 
rotation, for the A sublattice, is represented by the operator Ai. This is combined with the generators of rotations 
about the other two axes, and approximated over a finite interval in the following symmetric decomposition |14|: 



s(At) » A 3 ^At)A 2 ^At)A 1 (At)A 2 ^At)A 3 ^At)s 



(6) 



A similar approach is applied to the B-lattice updates. Each separate rotation is implemented with the finite rotation 
formula In the above equation, the rotation Ai associated with the largest eigenvalue Fji is centrally placed, 
and for brevity this is denoted the 32123 sequence; an empirical comparison with an alternative 12321 sequence is 
presented below. 

Ordering in this model is measured by the symmetric and traceless second-rank tensor 



1 



N 



Q = — . 
N ^ 

i=i 



2 S - 

2 J 



Note that this, like the hamiltonian of eqns and J2J, is invariant to all spin flips Sj — ► — Sj. The largest eigenvalue 
of Q is conventionally taken to be the order parameter Q, and the corresponding eigenvector n is the director. 

Comparison with Monte Carlo simulations is facilitated by evaluating the "configurational temperature" introduced 
by Rugh [l5j], and specifically derived for orientational degrees of freedom by Chialvo et al. [l^. The relevant 
expression, in our notation, is 



fc R T = 



i 



12 



(H) 



where Vj is the angular gradient, and V| the angular laplacian, with respect to the orientation of spin j. The 
second expression above is specific to the Lebwohl-Lasher potential, and is easily obtained from eqn J5J . In the results 
reported below, Boltzmann's constant fee, and the strength parameter J, are chosen to be unity. 

In Fig. H simulation averages generated by the spin dynamics algorithm are compared with those produced by 
conventional Monte Carlo. A system size of 8 x 8 x 8 spins is employed, which is sufficient to show interesting 
behaviour in the phase transition region, while being economical. It is worth emphasizing that the aim is not to 
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FIG. 1: Energy per spin, and nematic order parameter, as functions of temperature. Lines: Monte Carlo simulations. Circles: 
spin dynamics. Statistical errors are smaller than the plotting symbols. Filled symbols indicate state points studied in more 
detail below. 



properly characterize the transition, for which much larger systems are required 0, 0, 0. Spin dynamics runs of 
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20000 steps, each of length At = 0.01, starting from equilibrated Monte Carlo configurations, were carried out. Good 
agreement, within the statistical errors, is obtained with the Monte Carlo curves. This suggests that the sampling 
of the constant-energy hypersurface by spin dynamics generates satisfactory averages (but see below) . Three state 
points, at T ~ 1.5, 1.1,0.7, with order parameters Q w 0.0,0.5,0.8 representative of the disordered, transitional, and 
ordered states, respectively, were selected for further illustration. 

The accuracy of the algorithm, as measured by the root-mean-squared fluctuation of the energy 

AE Ima = v/(i? 2 ) - (£) 2 

is presented in Fig. |2J AE Ims cx At 2 over a wide range of At. For the largest timesteps there is some drift in the 
energy (which is removed in the calculation of AE ims ). Within each sublattice rotation step, there are several choices 
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FIG. 2: Root-mean-squared energy fluctuations, plotted against timestep, on logarithmic scales. Circles: T ~ 0.7, Q fa 0.8. 
Squares: T ~ 1.1, Q w 0.5 (displaced down by a factor of 10 for clarity). Diamonds: T fa 1.5, Q ~ (displaced down by a 
factor of 100 for clarity). Open symbols correspond to the rotation sequence 32123; filled symbols to the sequence 12321 (see 
text). The dashed line has gradient 2 for reference. 

for the sequence of axes about which to carry out the rotations. The figure shows that the sequence 12321 (in which 
the axis corresponding to the smallest eigenvalue of Fj is placed centrally) is worse than the sequence 32123 of eqn © 
(in which the largest eigenvalue is central) by a factor 2-3. Also the low-temperature ordered state point exhibits 
worse energy conservation than the high-temperature, disordered state point, reflecting the larger torques resulting 
from the local field. 

There are some caveats associated with spin dynamics applied to this model. Firstly, as noted, the total magneti- 
zation is conserved exactly by the dynamics, and asymptotically by the algorithm as At — > 0. This is an artificial, 
physically meaningless, conservation law for the present model. The hamiltonian, and the order parameter defined 
above, are invariant under all spin flips Sj — ► —Sj, reflecting the symmetry of the nematic phase. However, the dy- 
namics is not. If a spin is flipped, it will begin to rotate in the physically opposite direction. Macroscopic consequences 
come from this: in the ordered phase, the director rotates systematically about the fixed overall magnetization vector, 
if S = Yl j Sj is nonzero. This regular precession is superimposed on a general tendency of the director to become 
aligned with the magnetization vector: this happens slowly in the ordered phase, and more rapidly in the transition 
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region. The rate of precession is shown in Fig. El f° r both the state points, and for a range of tilt angles of the director 
relative to the magnetisation, at a range of net magnetisations obtained simply by flipping spins in an equilibrated 
Monte Carlo configuration. The precession rate is essentially independent of the tilt angle, and is proportional to both 




magnetisation per spin 

FIG. 3: Director rotation as a function of net magnetization per spin. Filled symbols: T ~ 0.7, Q ~ 0.8. Open symbols: 
T ~ 1.1, Q » 0.5. The director is inclined with respect to the magnetization vector by 30° (circles), 45° (squares), 60° 
(diamonds) and 90° (triangles). 

the order parameter Q, and to the net magnetization per spin, as would be expected by considering the interaction 
between a typical spin and the local field. 

This effect is a potential pitfall in simulations of these systems by spin dynamics. In a typical configuration of N 
spins, the net magnetization per spin will statistically be of order l/y/~N: this will produce a long-lived slow rotation 
of the director which, if uncorrected, will dominate measured dynamical properties. If the magnetization is, for any 
reason, substantial, the fast director precession distorts the measured simulation averages, such as configurational 
temperature. However, the solution to this problem seems relatively straighforward. Since the statistical properties 
of the Lebwohl-Lasher model are invariant to spin flips, it should be satisfactory to fix the system magnetization at 
the minimum possible magnitude, by flipping spins, at the start of a simulation run. In addition, once the director is 
aligned with the magnetization, the effects seem to be minimal. 

In this paper, an algorithm for simulating the Lebwohl-Lasher model by spin dynamics has been presented. It 
remains to be seen whether this approach will lead to the determination of interesting "dynamical" properties of this 
model, and related models, and possibly to accelerated sampling of the transition region. 
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